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Population structure induced by both spatial embedding and more general networks of interaction, such 

cn ; 

■ as model social networks, have been shown to have a fundamental effect on the dynamics and outcome of 

evolutionary games. These effects have, however, proved to be sensitive to the details of the underlying 

in ' 

topology and dynamics. Here we introduce a minimal population structure that is described by two distinct 
hierarchical levels of interaction, similar to the structured metapopulation concept of ecology and island 
Qh ■ models in population genetics. We believe this model is able to identify effects of spatial structure that 
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q-( the scope of our approach, we expect that those we are able to reproduce should be generally applicable 

to a wide range of models. We derive the dynamics governing the evolution of a system starting from 



do not depend on the details of the topology. While effects depending on such details clearly lie outside 



fundamental individual level stochastic processes through two successive meanfield approximations. In 
our model of population structure the topology of interactions is described by only two parameters: the 
effective population size at the local scale and the relative strength of local dynamics to global mixing. We 
demonstrate, for example, the existence of a continuous transition leading to the dominance of cooperation 
in populations with hierarchical levels of unstructured mixing as the benefit to cost ratio becomes smaller 
then the local population size. Applying our model of spatial structure to the repeated prisoner's dilemma 
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we uncover a novel and counterintuitive mechanism by which the constant influx of defectors sustains 
cooperation. Further exploring the phase space of the repeated prisoner's dilemma and also of the "rock- 
paper-scissor" game we find indications of rich structure and are able to reproduce several effects observed 
in other models with explicit spatial embedding, such as the maintenance of biodiversity and the emergence 
of global oscillations. 
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I. INTRODUCTION 



The dynamics of Darwinian evolution is intrinsically frequency dependent, the fitness of in- 
dividuals is tightly coupled to the type and number of competitors. Evolutionary dynamics acts, 
however, on populations, not individuals and as a consequence depends on not only population 
composition, but also population size and structure. Evolutionary game theory came about as the 
result of the realization that frequency dependent fitness introduces strategic aspects to evolution 
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IUI,I2L13D- More recently the investigation of the evolutionary dynamics of structured populations, 
where individuals only compete with some subset of the population, e.g. their neighbors in space 
or more generally in some graph J4, Q], has lead to the recognition that the success of different 
strategies can be greatly influenced by the topology of interactions within the population. Funda- 
mental differences were found - compared to well-mixed populations, where individuals interact 
with randomly chosen partners - in models that describe the evolution of cooperation (variants of 
the prisoner's dilemma game [4, (j, 7, 8, 2I) or deal with the maintenance of biodiversity in the 
context of competitive cycles (variants of the rock-paper-scissors game |3, Jl3, ill[3,[l3, 14]). 

In order to investigate the coevolutionary dynamics of games on structured populations the 
full set of connections between a potentially very large number of individuals must be specified. 
This is only possible by reducing the number of degrees of freedom considered, either through 



postulating a highly symmetrical (such as lattices 



i y, y, 0, h is 



random connection structure (such as some random graph ensemble [21 



9L |20fl) or fundamentally 



220). The question of how 



one goes about the task of reducing the number of degrees of freedom - of choosing the relevant 
parameters to describe the population structure constrained to which individuals undergo evolution 
- is not trivial. Both the explicit spatial as well as the random graph ensemble approach have clear 
precedents in condensed matter physics and network theory, respectively. It is not, however, clear 
which - if either - approach best describes natural populations of cyclically competing species or 
societies composed of individuals playing the prisoner's dilemma game. 

As an example let us consider colicin producing bacteria, that play the so called "rock-paper- 
scissors" (RPS) game (for details see below). This system has recently been the subject of two 
experimental studies aimed at demonstrating the role of structured populations in the maintenance 
of diversity. In the first study Jlp[ [lli bacteria were cultured in vitro in Petri dishes, effectively re- 
stricting competition between bacteria to neighbors on the (2D) Petri dish surface (Fig.l top left), 
while in the second experiment [12] in vivo bacterial colonies were established in co-caged mice 
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and their development was subsequently followed. In the case of the first experiment the analogy 
with explicit 2D spatial embedding (present by construction) is clear (Fig.l bottom left). The pop- 
ulation structure of the second experiment is, however, clearly different. The bacteria in individual 
mice can be readily considered as locally well-mixed populations, the coevolutionary dynamics 
of which reduces in the standard meanfield limit to a system of non-linear differential equations 
(the adjusted replicator equations I124ll ). As the experiments show, however, migration of bacteria 
between mice may also occur - resulting in the observed cyclic presence of the three strains in 
individuals. There are two distinct scales of mixing present in the system. Bacteria within each 
mice compete with each other forming local populations - an unstructured neighborhood com- 
posed of individual bacteria, while also being exposed to migrants from mice with whom they 
share the cage, together forming a global population - an unstructured neighborhood composed of 
individual local populations (Fig.l top and bottom right). This setup is referred to in the ecology 
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literature - albeit in significantly different contexts - as a "structured metapopulation" H26L 
where structured here refers to the detailed consideration of the population dynamics of the indi- 
vidual populations (often called "patches") comprising the metapopulation and is also related to 
the finite island models of population genetics |28|]. 

The above example of co-caged mice is not unique, we may readily think of other ecological or 
sociological examples where an approximation with hierarchical scales of mixing with no internal 
structure can be relevant (such as human societies with two distinct scales of mixing present, the 
first within individual nations the between them at an international level). We have, also, recently 
used a similar approach to construct a model of genetic exchange among bacteria of the same 
species (the bacterial equivalent of sex) with which we were able to take into account the effects of 
spatial and temporal fluctuations in a manner that can explain the benefit of such genetic exchange 
at the level of the individual 113 ill . 

In this paper we construct a hierarchical meanfield theory where the two distinct (i.e. local and 
global) scales of mixing are each taken into account in terms of two separate meanfield approx- 
imations and fluctuations resulting from finite population size on the local scale of mixing are 
also considered. We subsequently explore the similarities and differences between this and other 
models of structured populations in the case of the "rock-paper-scissors" and prisoner's dilemma 
games. Through these examples we suggest that our approach allows the separation of the effects 
of structured populations on coevolutionary dynamics into effects which are highly sensitive to 
and dependent on the details of the topology and those which only require the minimal structure 
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Petri dish experiment 



Co-caged mice experiment 





Explicit spatial embeding Two distinct scales of mixing 





FIG. 1: (Color online) In the colicin version of the RPS game, strains that produce colicins (red/dark 
grey) kill sensitive (green/light grey) strains, that outcompete resistant (blue/black) strains,that outcom- 



pete colicin producing strains (toxin production involves bacterial suicide). Experiments II 1 OTI show that 
colicin-producing strains cannot coexist with sensitive or resistant strains in a well-mixed culture, yet all 
three phenotypes are recovered in natural populations. Two recent experiments have examined the role of 



population structure in the maintenance of diversity among colicin-producing bacteria. In the first H 1QTI in 



vitro colonies were established on an agar substrate in Petri dishes, a setup which effectively limits compe 



tition to neighbors on the petri dish in analogy with explicit spatial embedding in 2D. In the second 01211 in 
vivo colonies were established in the intestines of co-caged mice, a setup which has two distinct scales of 
mixing, with no explicit structure on either scale. 



present in our approximation and can consequently (in terms of sensitivity to the details of the 
topology) be considered more robust. 
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II. HIERARCHICAL MEANFIELD THEORY FOR TWO DISTINCT SCALES 



Let us consider an evolutionary game between d types (strategies) described by the d x d payoff 
matrix A with elements a k j. Assuming finite and constant population size, natural selection can 



be described at the level of the individual by the so called the Moran process [30], during which at 
each time step an individual is selected randomly from the population to be replaced (death) by the 
offspring of an individual that is chosen proportional to its fitness to reproduce (birth). This models 
a population in equilibrium, where the time scale of the population dynamics is set by the rate at 
which "vacancies" become available in the population. The fitness of each individual depends on 
the payoff received from playing the game described by A with competitors (an individual of type 
k receiving a payoff a k j when playing with an individual of type j). In well-mixed populations, 
individuals can be considered to come into contact (compete) with equal probability with any 
member of the population excluding themselves - this allows one to calculate the fitness of an 
individual of type k in a meanfield manner, yielding 

^ a kj {rij - S kj ) 

7Tfc = 7T base + N _ 1 , (1) 

where n k is the number of individuals of type k in the population, J2k=i n k = N is the size of 
the population, 7r basc is some baseline fitness and the Kronecker delta symbol 5 k j is equal to unity 
if k = j and is zero otherwise. From this we may calculate the transition probabilities of our 
stochastic process, i.e., the probability of an individual of type i being replaced by an offspring of 
an individual of type k is given by 

_ m ir k n k 

lk N itN K ) 

where n = J2k=i n k n k/N. The state of any population is completely described by the frequency 
of the different strategies x k = n k /N. Due to the normalization J2 k=1 x k = 1, the values of x k 
are restricted to the unit simplex Sj yfl. For d = 2 this is the interval [0, 1], S 3 is the triangle with 



vertices {(1, 0, 0), (0, 1, 0), (0, 0, 1)} while SaJs a tetrahedron etc. 



As Traulsen et al. have recently shown 11241 12511 for sufficiently large, but finite populations 



the above stochastic process can be well approximated by a set of stochastic differential equa- 
tions combining deterministic dynamics and diffusion (population drift) referred to as Langevin 
dynamics: 

(2-1 

x k = a k (x) + y^Cfcj (x)£j (t) , (3) 

3=1 



where the effective deterministic terms afc(x) are given by 

a fc (x) = 2_^( T jk - T kj) = Xk -r-, , (4) 

Cfcj(x) are effective diffusion terms, that can also be expressed in terms of the transition probabil- 
ities as described in [fl], and £j are delta correlated (6fc(06j(O) = $kj$(t — f) Gaussian white 
noise terms. As N — > oo the diffusion term tends to zero as 1/VN and we are left with the 
modified replicator equation. 

In the context of our hierarchical mixing model the topology of connections can be described 
by two parameters, the populations size at the local scale of mixing N, and a second parameter 
/i, which tunes the strength of global mixing relative to the local dynamics. We take into account 
the second (global) scale of mixing - mixing among local populations - by introducing a modified 
version of the Moran process. In the modified process a random individual is replaced at each 
time step either with the offspring of an individual from the same population (local reproduction) 
or with an individual from the global population (global mixing). This is equivalent to considering 
the global population to be well-mixed at the scale of local populations. 

Let us consider a global population that is composed of M local populations of size N. In each 
local population vacancies become available that local reproduction and global mixing compete to 
fill. In any local population I the probability of an individual of some type k filling a new vacancy 
due to local reproduction must be proportional to the number of individuals of type k multiplied by 
their fitness i.e. Tv l k n l k , where we consider tx\ to be determined only by interactions with individuals 
in the same local population according to equation (OQ). To describe the tendency of individuals of 
some type k in local population I to contribute to global mixing we introduce the parameters a l k . 
The choice of appropriate a l k depends on the details of the global mixing mechanism, for systems 
where only the offspring of individuals mix globally it is proportional to the fitness of a given type, 
while for mechanisms such as physical mixing, by e.g. wind or ocean currents, it may be identical 
for each type. Irrespective of the details, however, the probability of an individual of some type 
k filling in a new vacancy due to global mixing should be proportional to the global average of 
the number of individuals of type k multiplied by their mixing tendency, which we denoted as 
(&kn k ) = YlaLi a k n k/M, and the strength of global mixing /i. These consideration lead to the 
new transition probabilities: 
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4fc ^ VELiWW + ^*«*» / * V JV(tH + /i(a)) ; ' 
where tx 1 = Y? k=1 n l k n l k /N and (a) = £jj =1 (<7 fc ra fc ) / 'TV '. 

We have found that the results presented below are qualitatively the same for both the fitness 
dependent choice of a\ = ir l k and the fitness independent choice of a l k — 1. Therefore, in the fol- 
lowing we restrict ourselves to the somewhat simpler fitness independent choice of oj. = 1, which 
can be considered to correspond to some form of physical mixing mechanism. The transition 
probabilities © then reduce to: 

fi = r A( ^ ^Wk , » ( n k) \ (6) 

ik N \n l + ii tt 1 N n l + fi N J K ' 

We can see that after a vacancy appears either local reproduction occurs, with probability tx 1 / [tx 1 + 
/i), or global mixing, with probability ji/{jx l + ji). From © we may derive the Langevin equation 
describing the coevolutionary dynamics of population I from the 

d-l 

x[ = a fc (x ; , (x» + ^c fci (x J , (x))£j(t), (7) 
with the modified deterministic terms given by 

A (v i _ 4KC^) - tt(x')) + //((a*,) - 4) 

, \X/J — , , {O) 

7T(X ( ) + /i 

where the vector (x) = YldLi^/M w i m components (x k ) = J2bLi x i/M describes the fre- 
quencies of the individual types in the global population and the diffusion terms c(x l , (x)) can be 
expressed in terms of the modified transition probabilities T\ k as above. 

Equations © describe the coevolutionary dynamics of the global population through the cou- 
pled evolution of the {x 1 , . . . , x M } local populations. In the limit of a large number of local popu- 
lations (M — > oo) the distribution of the local populations over the space of population states (the 
simplex Sd) is described by a density function p(x) that is normalized over Sd, i.e., f s p(x) = 1. 
The time evolution of p(x) follows ad—1 dimensional advection-diffusion equation - the Fokker- 
Planck equation corresponding to eq. ©: 

p(x) = -V{a(x, (x))p(x)} + iv 2 {b(x, (x»p(x)} , (9) 

with the global averages (x k ) = J s x k p(x) coupled back in a self-consistent manner into the 
deterministic terms ctfc(x, (x)) and the diffusion matrix b k j(x, (x)) = YltZi Cfci( x > (x))cy(x, (x)). 
For large local populations (N — > oo) the diffusion term vanishes as 1/JV. 
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The above advection-diffusion equation © presents an intuitive picture of the coevolutionary 
dynamics of the population at a global scale. We can see that local populations each attempt to 
follow the trajectories corresponding to the deterministic replicator dynamics, while under the 
influence of two additional opposing forces: (i) global mixing, which attempts to synchronize 
local dynamics and (ii) diffusion resulting from finite population size effects, which attempts to 
smear them out over the simplex. The strength of these forces are tuned by two parameters /x and 
N, respectively. 

If, further, the effects of synchronization are irrelevant, as for example in the case of populations 
where selection is externally driven by independent environmental fluctuations, we may replace 
the global population average with the time average of any single population. This is the approach 



we used in our study of genetic mixing in bacteria 113 111 . 

During our numerical investigations we found solving the advection-diffusion equation © nu- 
merically challenging, particularly in the N — > oo limit. We resorted instead to solving the coupled 
Langevin equations © for large M = 10 4 — 10 5 to simulate the time evolution of p(x). 

IE. COOPERATION IN POPULATIONS WITH HIERARCHICAL LEVELS OF MIXING 

The evolution of cooperation is a fundamental problem in biology, as natural selection under 
most conditions favors individuals who defect. Despite of this, cooperation is widespread in na- 
ture. A cooperator is an individual who pays a cost c to provide another individual with some 
benefit b. A defector pays no cost and does not distribute any benefits. This implies the payoff 
matrix 

b — c —c 
b 

where b is the benefit derived from playing with a cooperator while c is the cost for cooperation. 
From the perspective of evolutionary game theory, which equates payoff with fitness, the apparent 
dominance of defection is simply the expression of the fact that natural selection a priori selects 
for fitness of individuals and not the fitness of groups. 

Defection dominates cooperation in any well-mixed population [3]. Population structure in- 



(10) 



duced by spatial structure [|4|, UM and more general networks of interactions [|2lL I22L |23|]) has, 



however, been found to facilitate the emergence and maintenance of cooperation. The mecha 



nism responsible, termed spatial, or more generally, network reciprocity[32] depends strongly on 
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FIG. 2: a In an infinitely large well-mixed population evolutionary dynamics is deterministic and leads 
to the extinction of cooperators as average fitness monotonically declines. The only stable fixed point 
corresponds to the point where the fraction of cooperators is zero (x = 0). To understand qualitatively the 
mechanism favoring cooperation in hierarchically mixed populations let us consider some density of local 
populations (p(x, t)) that is symmetric around its mean at time t. Due to global mixing all local populations 
are being driven toward the global average. Due to the influx bias, however, populations with a lower than 
average number of cooperators will be driven stronger (faster) than those on the other side of the average. 
Examining the density of local populations at some time t + At, this results in a net movement of the global 
average toward a larger fraction of cooperators. This is, of course, opposed by local reproduction that favors 
an increase in the number of defectors. For the global average to keep moving toward a higher number of 
cooperators and eventually to keep balance with local reproduction bias a density of local population with 
finite width is needed over which the effect of the influx bias can exert itself. It is drift caused by local 
population size that maintains this finite width, and this is the reason that the b/c threshold above which 
cooperation dominates depends on local population size, b Stationary density of local populations p(x) for 
different values of b/c with N = 100, [i = 0.1. c Transition toward a global dominance of cooperation for 
/i = 10. (triangles), fi = 1 (crosses), n = 0.1 (squares), /x = 0.01 (circles) with N = 100. The critical value 
of b/c depends only weakly on [i changing by 20% over four orders of magnitude d Critical values of b/c 
as a function of N for different values of fi (notation as before). The dashed line corresponds to b/c = N. 
The critical b/c values were determined by numerically finding the inflection point of the transition curves. 
M = 10 3 was used throughout. 
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the details of local topology. In particular, it seems that lattice like connectivity structures where 



three- site clique percolation occurs II 1711 and more general interaction graphs where the degree of 



nodes k does not exceed the ratio of benefit to cost (i.e. k < b/c) I122ll are required for cooperation 
to be favored. 

Examining the effects of hierarchical mixing on the evolutionary dynamics of cooperation we 
found that a sharp, but continuous transition leads to the dominance of cooperation as the benefit 
to cost ratio becomes smaller then the local population size, i.e. b/c < N. If the benefit to cost 
ratio is larger then the local population size the global population is dominated by defectors. The 
mechanism leading to the dominance of cooperation arises due to the competition between local 
reproduction and global mixing. In local populations with lower average fitness - larger number of 
defectors - the influx of individuals from the global scale will be larger than in local populations 
with higher average fitness (cf. eq. (6) where the relative strength of the two terms on the left hand 
side depends on the sum of the average fitness of population I and ji). The crucial ingredient for 
cooperation to be successful is population drift introduced by finite local population size. It is bi- 
ased influx coupled with drift that can result in cooperation being favored in the global population 
(FigEl). 

IV. THE RPS GAME 

To explore the effects of hierarchical mixing in the context of games with three strategies we 
first turn to the case of the so called "rock-paper-scissors" (RPS) game. In the original popular 
version of the game two players are afforded the chance to simultaneously display either rock 
(fist), paper (flat hand) or scissors (two fingers). If player one displays a flat hand while player 
two displays a fist, player one wins as paper wraps rock. Similarly scissors cut paper, and rocks 
smashes scissors. Several examples of this game have been found in nature (e.g. among lizards 



|33|] ), but it is bacteria that have received the most experimental and theoretical attention. 

In ecology the often high diversity among microbial organisms in seemingly uniform envi- 
ronments, referred to as the "paradox of the plankton", has been difficult to understand. Several 



models 



?ased on spatially explicit game theoretical models have been proposed to explain this di- 



versity 



ityfllll, 
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14|1 . These models are all variants of the RPS game played by colicin producing 
bacteria. Colicins are antibiotics produced by some strains of Echerichia coli. In experiments (see 
Fig.l) typically three strains are used: colicin producing (C), sensitive (S) and resistant (R). The 
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FIG. 3: a Deterministic replicator dynamics (the N — > oo limit) of the symmetric RPS game consists of 
neutrally stable orbits along which the product of the strategy frequencies x^xpxs is conserved. If global 
mixing is present > 0) local populations deviate from these neutral orbits toward the global average 
(x) . Considering the simplest system with global mixing, that consisting of M = 2 local populations we 
see that in the presence of global mixing population x 1 and population x 2 move toward each other, respec- 
tively moving closer and further from the barycentre of the triangle until they become synchronized and 
subsequently pursue a common orbit. For deterministic local dynamics (N — > oo) such synchronization 
invariably occurs for any M if fj, > and typically converges to the barycentre of the simplex for suffi- 
ciently homogeneous initial conditions, b The deterministic replicator dynamics of the repeated PD game is 
markedly different from that of the RPS game in that the internal fixed point is unstable and in the absence 
of global mixing only ALLD survives. Again turning to the simplest scenario with M = 2 we see that 
if \i = any pair of populations x 1 and x 2 (gray and black lines) converge to the to the ALLD corner. 
As is increased above a critical value a second, stable configuration emerges: for a large subset of the 
possible initial conditions (all, but the left most x 2 ) we see that one of the populations (x 1 ) converges to 
ALLD , while the second (x 1 ) approaches a limit cycle. If fi is increased further, the above configuration 
ceases to be stable, the population which initially converges to ALLD (x 1 ) is subsequently "pulled out" 
by global mixing, following which the two populations synchronize and are finally absorbed together in 
ALLD. Simulations, however, show that synchronization may be avoided for M > 2 if fi is not too large. 

coevolutionary dynamics of the three strains can be cast in terms of an RPS game, C strains kill S 
strains, but are outcompeted, by R strains, because toxin production involves the suicide of bacte- 
ria. The cycle is closed by S strains that outcompete R strains, because resistance requires mutant 
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versions of certain membrane protein, which are less efficient than the wild type fl 1 ON - Despite 
the cyclic dynamics colicin-producing strains cannot coexist with sensitive or resistant strains in a 
well-mixed culture, yet all three phenotypes are recovered in natural populations. Local dispersal 
(modeled as explicit spatial embedding) has widely been credited with promoting the maintenance 
of diversity in this system 

In its most symmetric form the RPS game is described by the payoff matrix 



/ 
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and some 7Tb a sc > e- The dynamics of this game in an infinitely large well mixed population 
consists of neutral orbits along which the product x^xpXs is conserved. For any finite N, however, 
fluctuations lead to the inevitable extinction of all but one of the strategies 111511 . Spatial population 
structure can avert this reduction in diversity fl 1 CH. Il3ll through the emergence of a stable fixed 
point at the barycentre of the simplex . The effect of the gradual randomization of different lattice 
topologies (where a small number of edges are randomly rewired) on the dynamics of the game 
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has also been investigated. A Hopf bifurcation leading to global oscillations was observed [34, 
as the fraction of rewired links was increased above some critical value. 

Examining the dynamics of the symmetric RPS game in terms of our hierarchical meanfield 
approximation we observed that an internal fixed point emerged for N — > oo (Fig. 3 a). More 
importantly, diversity was also maintained for finite local population sizes if global mixing was 
present. Simulations of the time evolution of p(x) also revealed a Hopf bifurcation leading to 
the oscillation of the global average as fi was increased above a critical value /i c depending on 
N (Fig. 4a). These results show that previous results obtained from simulations of populations 
constrained to different lattice topologies can be considered universal in the sense that not only 
lattices, but any population structure that can be approximated by two distinct internally unstruc- 
tured scales of mixing are sufficient for their existence. In the context of the "paradox of the 
plankton" these results imply that aside of local dispersal (modeled as explicit spatial embedding) 
a minimal metapopulation structure (with local competition and global migration) can also facili- 
tate the maintenance of diversity in cyclic competition systems. 
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FIG. 4: (Color online) a In the case of the rock-paper-scissors game a 



observed for populations evolving on gradually randomized lattices H34 . 



Hopf bifurcation similar to that 



3511 leads to the emergence of 



global oscillations (the red line indicates the trajectory of (x)) if /i is larger than a critical value fJ, c {N) (see 



video S 1 [36]). The density p(x) is indicated with a blue color scale, b The ratio A of the area of the global 
limit cycle and the area of the simplex is plotted as a function of /i for three different values of N. For the 
repeated prisoner's dilemma game the combination of finite local population size and global mixing \i > 
can lead to a stationary solution (c) qualitatively similar to that observed for explicit spatial embedding e. 
This state is characterized by a stable global average (large dot), just as the lattice system (data not shown) 
and sustained local cycles of cooperation, defection and reciprocity, also similar to the lattice case where 
groups of ALLD (red, dark grey) individuals are chased by those playing TFT (blue, black), which are 
gradually outcompeted by ALLC (green, light grey), d As \i is decreased a discontinuous transition can be 
observed to the ALLD phase. The ratio I of populations on the internal cycle is plotted as a function of 
[i. The inset shows the transition for different values of N. f The same critical line in the fj,-N plane can 
be approached by increasing N with /i fixed. A large hysteresis can be observed as N is decreased below 
the critical value indicating the discontinuous nature of the transition. We numerically simulated the time 
evolution of p(x) by integrating the stochastic differential equation system denned by eq. ^} for large M 
(10 4 - 10 5 ) throughout. For the RPS game we used 7Tb a se = 1 and e = 0.5, while in the case of the repeated 



PD game we followed ref. Q38Q, setting T = 5, R = 3, P = 1, S = 0.1, m = 10 and c = 0.8. Lattice 
simulations (e) where performed on 1000 x 1000 square lattice with an asynchronous local Moran process 
between neighbors and periodic boundary conditions. 
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V. THE REPEATED PRISONER'S DILEMMA GAME 



In the general formulation of the prisoner's dilemma (PD) game, two players have the choice 
to cooperate or to defect. Both obtain some payoff R for mutual cooperation and some lower 
payoff P for mutual defection. If only one of the players defects, while the other cooperates, 
the defector receives the highest payoff T and the cooperator receives the lowest payoff S. That 
is T > R > P > S and defection dominates cooperation in any well-mixed population. New 
strategies become possible, however if the game is repeated, and players are allowed to chose 
whether to defect or cooperate based on the previous actions of the opponent. In the following 
we consider, similar to refs. D7I1 and 113 811 that recently examined the role of finite population 
size and mutation and finite population size, respectively in terms of the repeated PD game with 
three strategies: always defect (ALLD), always cooperate (ALLC), and tit-for-tat (TFT). TFT 
cooperates in the first move and then does whatever the opponent did in the previous move. TFT 
has been a world champion in the repeated prisoner's dilemma ever since Axelrod conducted his 
celebrated computer tournaments 07J], although it does have weaknesses and may be defeated by 
other more complex strategies [|39ll . 

Previous results indicate that if only the two pure strategies are present (players who either 
always defect or ones who always cooperate) explicit spatial embedding [4] and some sufficiently 
sparse interaction graphs 11221 14011 allow cooperation to survive and the behavior of populations is 



highly sensitive to the underlying topology of the embedding [17]. We have found that introducing 
global mixing into the PD game with only the two pure strategies present also allows cooperation 
to survive. The mechanism responsible for favoring cooperation in this case, however, depends 
on the details of the competition between local reproduction and global mixing. For more than 
two strategies these details are much less relevant and do not qualitatively influence the dynamics. 
We will, therefore, consider the delicate issues concerning the PD game with only the two pure 
strategies in a separate publication, and concentrate here on the repeated PD game with three 
strategies. 

To investigate the effect of global mixing on the repeated PD game with three possible strate- 



gies: ALLD, ALLC and TFT following Imhof et al. 



3811 we considered the payoff matrix: 
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deterministic local populations 



t 



(i) 



(ii) 



(iii) 



N 



effective single population 



(i) 




FIG. 5: (Color online) Phase space for the repeated prisoner's dilemma game on a population structure 



with two distinct scales (see video S2 j36tl ). Three different phases are possible depending on the values of 
fj, and N: (i) only ALLD survives (ii) an internal limit cycles is maintained by global mixing due to a large 
density of local populations around the ALLD corner (iii) a globally oscillating self maintaining limit cycle 
is formed. For extreme values of /x the global dynamics reduces to that of some well-mixed population 
where only ALLD survives: As fi becomes negligible (// <C ir^ for all k) we approach the limit of isolated 
local populations, while for /i > ^ we are left with a single synchronized population. Similarly for N = 2 
- the smallest system with competition - the system can be described as a single well mixed population for 
any n and ALLD again prevails. In the limit of deterministic local populations (N — ► oo) all three phases 
can be found depending on the value of fi. The density p(x) is indicated with the color scale. A figure 
illustrating the phase space of the repeated prisoner's dilemma game with fitness dependent global mixing 



is included in the supplementary material 13611 . 
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where the strategies are considered in the order ALLC, ALLD, TFT, m corresponds to the number 
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of rounds played and c to the complexity cost associated with conditional strategies (TFT). The 
dynamics of this game has a single unstable internal fixed point and the state where each member 
of the population plays ALLD is the only nontrivial stable equilibrium (Fig. 3b). 

Introducing global mixing, between local well-mixed populations, however, causes new sta- 
tionary states to emerge . Three phases can be identified: (i) ALLD wins (ii) large fraction of local 
populations in the ALLD corner maintains local cycles of cooperation defection and reciprocity 
through providing an influx of defectors that prevent TFT players from being outcompeted by 
ALLC playing individuals (iii) a self maintaining internal globally oscillating cycle emerges. The 
simplest scenario of two (M = 2) deterministic (N — > oo) local populations coupled by global 
mixing (fi > 0) already leads to the emergence of phase (ii) as demonstrated in Fig fib while phase 
(iii) only emerges for larger M. For larger M simulations show that in the limit of large local pop- 
ulations all global configurations with less than some maximum ratio of the populations 7 on the 
internal cycle are stable in phase (ii). A transition from phase (ii) to (i) happens as fi is decreased 
below a critical value /i"^ 1 and 7 approaches zero as 7 = (1 — /i"^ 1 ///) (data not shown). This can 
be understood if we considered that near the transition point a critical proportion C = fj,(l — 7) of 
ALLD individuals needs to arrive to stabilize local cycles of cooperation defection and reciprocity. 
At the critical point 7 = and ii = /i"^ 1 which implies C = ji 1 ^ 1 giving 7 = (1 — ll 1 ^ 1 / ' ll) 

Exploring the N — fx phase space (Fig. 5) we see that the transition from phase (i) to (ii) be- 
comes discontinuous for finite N (Fig.4d,e). Further, for any given value of N and fi the global 
configuration is described by a unique 7 due to the presence of diffusion. For appropriate values 
of the parameters the global average converges to a stationary value in phase (ii) similarly to case 
of explicit spatial embedding (Fig. 4c). 

For very small (p <C ix k for all k) and very large (p, <C 7r fc ) values of it, the global dynamics can 
be reduced to that of some well-mixed population where only ALLD persists (Fig 5.). For small N 
we again have an effective well-mixed population - the only limit were defectors do not dominate 
is N — > oo. In comparison with previous results of Imhof et al. we can see that evolutionary 
cycles of cooperation defection and reciprocity can be maintained not only by mutation, but also 
by population structures with hierarchical levels of mixing. 



16 



VI. DISCUSSION 



While it is, of course, clear that the reduction of any realistic population structure to a manage- 
able construction is always an approximation, it has not been clearly established what the relevant 
degrees of freedom are in terms of evolutionary dynamics. Meanfield approximations are a classic 
method of statistical and condensed matter physics and are routinely used to circumvent intractable 
combinatorial problems w 



/hich arise in many-body systems. Cluster- meanfield approximations of 
sufficient precision lll8L ll9ll have been developed that adequately describe the evolutionary dynam- 
ics of explicitly structured populations through systematically approximating the combinatorial 
complexity of the entire topology with that of small motif of appropriate symmetry. The effects of 
more minimal effective topologies have, however, not been investigated previously. In the above 
we have shown that straightforward hierarchical application of the meanfield approximation (the 
assumption of a well-mixed system) surprisingly unveils a new level of complexity. 

In the broader context of ecological and population genetics research on structured populations 
our model can be described as a metapopulation model. The term 'metapopulation' is, however, 
often used for any spatially structured population Ii27ll . and models thereof. More restrictive defi- 
nitions of the term are often implied in the context of ecology and population genetics literature. 

The foundations of the classic metapopulation concept where laid down by Levin's vision of 
a "metapopulation" as a population of ephemeral local populations prone to extinction. A classic 
metapopulation persists, like an ordinary population of mortal individuals, in a balance between 
'deaths' (local extinctions) and 'births' (establishment of new populations at unoccupied sites) 
Il27n . This classic framework is most wide spread in the ecology literature, a less often employed 
extension is the concept of a structured metapopulation where the state of the individual popula- 
tions is considered in more detail, this is more similar to our concept of hierarchical mixing, but 
differs in considering the possibility of local extinctions. 

The effects of finite population size and migration, which our model considers, has been of 
more central concern in the population genetics literature. The analog of Levin's classic metapop- 
ulation concept is often referred to as the 'finite-island' mod el Il28n the effective population genetic 
parameters describing which, have been explored in detail Ii29ll . The study of the population ge- 
netics of spatially subdivided populations in fact predates Levin, Wright having emphasised the 
capacity of drift in small populations to bring about genetic differentiation in the face of selection 
and/or migration several decades prior[|280. 
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Our hierarchical mixing model treats the coevolutionary dynamics of evolutionary games on 
structured populations in a manner similar to the most simple population genetic models of spa- 
tially subdivided populations, focusing on the parallel effects of selection, drift and migration. It 
goes beyond these models both in considering the effects of frequency dependent selection (and 
the strategic aspects of the evolutionary dynamics this implies) and in using a self-consistent ap- 
proach to describe the global state of the subdivided population. Also, in order to maintain a 
connection with previous work on the effects of spatial structure on evolutionary games, which 
rely on Nowak's concept of spatial games flj], with individuals restricted to interact, and hence 
compete, only with neighbours as defined by some topology of interaction, we develop our model 
from the level of the individual by introducing a modified version of the Moran process - and not 
by extending the Wright-Fisher process (which considers discrete generations and binomial sam- 
pling to account for finite population size). The effective population structure described by our 
hierarchical mixing model can be thought of as a population of individuals, interactions among 
which are specified by the edges of a hierarchically organized random graph. The fundamental 
difference in our picture is that the edges of this graph of interactions are not considered to be 
fixed, but are instead in a constant state of change, being present with a different probability be- 
tween pairs of individuals who share the same local population and between pairs of individuals 
who do not (Fig.l.). We consider annealed randomness, which in contrast to the usual quenched 
picture of fixed edges is insensitive to the details of topology. Our approach we believe best facil- 
itates the exploration of the effects of changing the relative strengths of drift and migration in the 
context of evolutionary games on structured populations. 

Examining the effects of hierarchical mixing in the context of the evolution of robustness we 
demonstrated that biased influx coupled with drift can result in cooperation being favored, pro- 
vided the ratio of benefit to cost exceeds the local population size. This result bears striking 
resemblance to that of Ohtsuki et al. [12211 . who were able to calculate the fixation probability of a 
randomly placed mutant for any two-person, two-strategy game on a regular graph and found that 
cooperation is favored provided the ratio of benefit to cost exceeds the degree of the graph. Our 
results demonstrate that this rule extends to the minimal spatial structure induced by hierarchical 
levels of mixing. 

Applying our model of spatial structure to the repeated prisoners dilemma revealed that a con- 
stant influx of defectors can help to stabilize cycles of cooperation, defection, and reciprocity 
through preventing the emergence of an intermittent period of ALLC domination in the popula- 
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tion, which would present a situation that "leaves the door wide open" to domination by defectors. 
While previous work has been done on the effects of "forcing" cooperation [41] the idea that an 
influx of defectors can in fact stabilize the role of reciprocity in promoting cooperation has not 
been proposed previously. It seems highly unlikely that this mechanism can be explained in terms 
of kin or multilevel (group) selection, the similarities between which in structured populations 
have recently been the subject of intensive debate (see e.g. [42] and J43] or [44] and [45]). Kin 
selection can operate whenever interactions occurring among individuals who share a more recent 
common ancestor than individuals sampled randomly from the whole population H5\\ are relevant. 
In our case it is the interaction between defectors, arriving from the global scale, and TFT players 
present at the local scale that is important, and not the interaction between individuals in the local 
population, who may be thought of as sharing a recent common ancestor due to local dispersal. 
Also, while the concept of multilevel selection presents a promising framework for the study of 
evolution of cooperation, it must nonetheless be possible to derive it from "first principles" - just 
as kin selection can be cast as an emergent effect of local dispersal. 

While there has been considerable work on studying the evolutionary games on graphs and 
highly symmetric spatial structures very little attention has been paid to the effects of more min- 
imal effective population structures, despite their widespread application in ecology and popula- 
tion genetics, fields from which evolutionary game theory was born and must ultimately reconnect 
with. We believe that the minimal population structure that such a hierarchical meanfield theory 
describes is potentially more relevant in a wide range of natural systems, than more subtle setups 
with a delicate dependence on the details and symmetries of the topology. We showed through 
two examples that such structure is sufficient for the emergence of some phenomena previously 
only observed for explicit spatial embedding, demonstrating the potential of our model to identify 
robust effects of population structure on the dynamics of evolutionary games that do not depend 
on the details of the underlying topology. The practical advantage of our approach, lies in its 
ability to readily determine whether or not some feature of a structured population depends on the 
topological details of local interactions. 

Recent simulation result concerning the dynamics of public goods games on different popula- 

un 

tion structures fl9J,|46|] and experiments where global mixing in an RPS like bacteria-phage system 
lead to the emergence of a "Tragedy of the commons" scenario Ii47l1 should all be amicable to 
analysis in terms of our method. 
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VIII. APPENDIX 

Our approach readily generalizes for an arbitrary number of hierarchical mixing levels. For 
three levels of mixing we may consider the global population to be comprised of M. subpopula- 
tions each of which is in turn subdivided into M local populations. With m E {1, ■ ■ • , Ai} running 
over subpopulations and / G {1, • • • , M} over local populations the transition probabilities can be 
written as: 

fm i = nf( + j^W)v + ^ 2) ((™*? l ')i')m> \ 13 

^ N VLU{K l < + +^H(W*f l ')l')m>)J ' 

where primed indices indicate the scale of mixing over which the average is taken, ^ describes 
the strength of mixing, and the ^a™ 1 the tendencies of mixing among local populations within a 
subpopulation, while fi^ describes the strength of mixing, and the ^a™ 1 the tendencies of mixing 
among subpopulations in the global population. 
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